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We present a self-contained formalism for calculating the background solution, the lin¬ 
earized solutions, and a class of generalized Frobenius solutions to a system of scale 
invariant differential equations. 

We first cast the scale invariant model into its equidimensional and autonomous 
forms, find its fixed points, and then obtain power-law background solutions. After 
linearizing about these fixed points, we find a second linearized solution, which provides 
a distinct collection of power laws characterizing the deviations from the fixed point. 
We prove that generically there will be a region surrounding the fixed point in which 
the complete general solution can be represented as a generalized Frobenius-like power 
series with exponents that are integer multiples of the exponents arising in the linearized 
problem. This Frobenius-like series can be viewed as a variant of Liapunov’s expansion 
theorem. 

As specific examples we apply these ideas to Newtonian and relativistic isothermal 
stars and demonstrate (both numerically and analytically) that the solution exhibits 
oscillatory power-law behaviour as the star approaches the point of collapse. These series 
solutions extend classical results; as exemplified for instance by the work of Lane, Emden, 
and Chandrasekhar in the Newtonian case, and that of Harrison, Thorne, Wakano, and 
Wheeler in the relativistic case. 

We also indicate how to extend these ideas to situations where fixed points may not 
exist — either due to “monotone” flow or due to the presence of limit cycles. Monotone 
flow generically leads to logarithmic deviations from scaling, while limit cycles generally 
lead to discrete self-similar solutions. 

Keywords : power law; Frobenius; Liapunov; scale invariant; stellar structure 


1. Introduction 

The presence of power-law behaviour in nature is such an extremely common phe¬ 
nomenon that considerable lore has now grown up concerning its genesis. One of 
the most common situations in which it occurs is in the presence of scale-invariant 
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systems. Such behaviour occurs, for instance, in any sort of thermodynamic system 
undergoing a second-order phase transition, and the behaviour of physical quanti¬ 
ties (such as susceptibilities) in terms of the distance from criticality is typically 
given by a power-law: 


X °c 


T — ^critical \ ' 
^critical J 


(1) 


Second-order phase transitions are extremely well-described by statistical field the¬ 
ories. The associated technical machinery (the renormalization group), is now so 
well developed that it is sometimes difficult to remember that second-order phase 
transitions are not the only route to power-law behaviour. The onset of power-law 
behaviour actually occurs at a much more primitive level, and can be analyzed 
directly in terms of the underlying differential equations. 

Scale invariance is commonly thought to be synonymous with power-law be¬ 
haviour. This is not quite correct: Scale invariance is not enough to guarantee 
power-law behaviour, while power-law behaviour implies scale invariance only in 
the trivial sense that it is always possible to write down some scale-invariant dif¬ 
ferential equation. (But such a DE does not necessarily have to describe a physical 
theory.) In the present paper, we examine scale invariant theories and find one 
direct route to power-law behaviour. 

This formalism consists of transforming a scale invariant theory into its associ¬ 
ated autonomous forma and then determining the location of its fixed points. These 
fixed points will lead inexorably to a background solution that exhibits power-law 
behaviour. By expanding the system about its fixed point, a linearized autonomous 
differential equation emerges, whose solution will exhibit different power-law be¬ 
haviour. This naturally leads to two independent power laws: one for the back¬ 
ground fixed point and a second one for the linearized deviations from the fixed 
point. Combining both results, our formalism indicates that the complete general 
solution to the system will, over some limited region, have a Frobenius-like power- 
series expansion. The usual Frobenius series is one of the form 


OO 

y{x) =x?-Y j ai x l . 

i—0 


( 2 ) 


with i a natural number, and is suitable for exploring linear ODEs in the region 
of regular singular points. In the present context, involving nonlinear DEs, the 
generalized Frobenius-like series we encounter are generically of the form 


y{x) = x p ■ 




( 3 ) 


»=0 


The exponent in the prefactor of this Frobenius-like series (p) is governed by the 
exponent associated with the fixed-point power-law background behaviour, while 
the power series is given in terms of exponents (*A) that are integer multiples of 
those arising in the linearized problem (A). This type of expansion is essentially a 
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modification of that given by Liapunov’s expansion theorem for autonomous deJ, 
wherein we have now back-tracked to give the expansion in terms of the original 
scale-invariant system. 

As specific examples we apply this formalism to both Newtonian and relativistic 
isothermal stars, where the relevant differential equations are equivalent to the ex¬ 
ponential [isothermal] Lane-Emden equation and its relativistic generalization. The 
resultant Frobenius-like power series generalize the approximate solutions found, for 
instance, by Chandrasekhara and by Harrison et. al .□, and provides significantly 
greater analytic information concerning their structure. This formalism is also ap¬ 
plied to a more general star, where the equation of state need not be linear except 
in the core. Both Newtonian and relativistic stars are verified to undergo damped 
oscillations in the total mass, pressure and radius as the central density goes to 
infinity. These oscillations can accurately be described by a Frobenius-like power 
series with complex exponents. The relevant exponents (though not the numeri¬ 
cal values of the coefficients) drop out naturally from the linearized problem. This 
power law solution matches the behaviour observed in numerical calculations to 
great accuracy. 


2. Why power laws are ubiquitous 

In this section we lay down the basic foundations of the formalism. First, we define 
“scale invariant”, “equidimensional”, and “autonomous”, and present the appro¬ 
priate substitutions to go from one form to another. A single nth-order differential 
equation is chosen as a model because of its simplicity, and then the formalism is ex¬ 
tended to a system of I coupled nth-order differential equations. After determining 
the existence of fixed points in the equation, the background solution is obtained. 
The differential equation is then expanded about its fixed points to obtain a second 
power-series solution near the fixed point (linearized solution). Combining both re¬ 
sults, a Frobenius-like series solution is obtained in a neighborhood close enough to 
the fixed point. 


2.1. Fundamental definitions 

Let us define scale invariant differential equations in the usual manner, as those 
which remain unchanged under the substitutions x —> ax and y{x) —> a p y(ax ), 
where a p is a fixed but arbitrary parameter.!!! Equations that remain unchanged 
under the substitutions x —> bx and y(x) —> y(x), where b is another fixed but 
arbitrary parameteiu, will be defined as equidimensional-in-x. Thus equidimensional 
equations can always be viewed as scale invariant equations corresponding to the 
particular value p = 0 for the exponent. Finally, autonomous equations will be 
defined as those equations which remain unchanged under the substitutions x —> 
x + c and y(x) —> y(x), where c is a third fixed but arbitrary parameter!!]. This last 
condition is that of translation invariance in the independent variable x. 
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It is a standard result that it is always possible to transform scale invariant 
equations into equidimensional-in-x equations through the following substitution^ 

y(x) = x p ■ w(x). (4) 

It is also a standard result that we can transform equidimensional-in-a; equations 
into autonomous equations by substituting^ 

x = e t z(t) = w( x), (5) 

where t is the new variable. 


2.2. From scale invariance to autonomy 

Let us first consider a nth-order scale invariant differential equation in a single 
dependent variable, 

F ^x, y(x), y\ x), y"{ x ),( 6 ) 


where the primes stand for differentiation with respect to x. We will write this as 
F(x, y(x)) = 0, but the presence of an appropriate number of derivatives should be 
inferred. Now assuming this equation to be scale invariant, let’s apply equation (|]) 
to make equation (JgJ) equidimensional-in-x. The derivatives of y(x) become 


d n y(x) 
dx n 



[x p w(x)] 


y _ 

*—•' i! (n — 


n\ p ! 


y.p-0 


3=0 


j'. (n - j)\ (p - j)\ 


(x). 


(7) 


Note that the jth derivative of y{x) involves only the jt\i and lower-order derivatives 
of w(x). Rewriting F , we can define a new function F as follows, 


F ( x, w(x), w/(x),... 


d n u;(x)^ _ 
dx” 


F j x, x p w(x),px p : w(x) + x p w'(x ),..., y 

3=0 


n\p\ 


j\{n- j)\{jp- j)\ 


( 8 ) 

xP-Ju,(»-j)(x) ] , 


since F depends only on x, w(x) and its derivatives. The new differential equation 
F(x,w(x)) = 0 is now equidimensional-in-x. 

Let’s now apply the substitution in equation (|j) to transform it into an au¬ 
tonomous equation. The derivatives of z are 


d" 

dx" 



(9) 


As before, note here that the product in equation (jt]) involves an nth order term 
plus lower-order derivatives (and no higher-order derivatives). Rewriting F we can 














February 7, 2008 2:12 WSPC/Guidelines-IJMPA power-law 


Power laws, scale invariance, and generalized Frobenius series... 5 

define a new function F via 

F^z{t),z{t),...,^z{t)^ f(t) = 



where the dots denote differentiation with respect to t and /(f) is some arbitrary 
function of f. The function F(z(t)) is now translation invariant, since the substi¬ 
tution t —> t + b as z(t) — > z{t + b) leaves the equation F(x,w(x)) = 0 (based on 
setting the RHS = zero) unchanged. The only way this can be true is if we can 
factorize F(x,w(x)) into an autonomous F times some /(f), yielding the LHS. 

Our original scale invariant equation F(x, y{x)) = 0, that had become equivalent 
to the equidimensional-in-x equation F(x,w(x)) = 0, has now become equivalent 
to the autonomous equation F(z(t)) = 0. We can of course also work backward 
and prove that for any autonomous-in-t equation for the dependent variable z(t) 
the substitution x = e* (that is, d/dt = x A/Ax) produces an equation that is 
equidimensional-in-x with z(t) = z(lnx) = w(x). Furthermore for any value of 
p the substitution w(x) = x~ p y(p) results in a differential equation that is scale 
invariant. However, this scale invariant equation need not be equal to the physical 
differential equation that describes the system we are concerned with. Without 
some additional information beyond the mere existence of a power law, there is 
simply not much one can say about the physics. 


2.3. Scale-invariant ordinary differential equations of nth-order 
2.3.1. Background solution 

In order to obtain the background solution of the ODE, we must look at the fixed 
points. Let’s assume that F has a fixed point z * such that all derivatives of z(t) 
vanish at that fixed point, i.e., 

f(, W = z.,^=o)=0, (ID 

where i € [1, n] denotes any positive integer in that range. Hence, a solution to the 
differential equation F(z(t)) = 0 is z(t) = z *. Tracing back all the substitutions we 
made, we obtain 

y(x) = x p z* which solves F(x,y(x)) = 0. (12) 


From this equation, the power law behaviour is obvious. However, not all scale 
invariant equations possess fixed points. A scale invariant differential equation, 
whose associated autonomous equation does not have any fixed points, will not 
possess power-law solutions. Refer to Appendix A for an explicit counter-example. 
In light of this, we can formulate the following theorem: 
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Theorem 1: Any first-order differential equation that is scale invariant and whose 
associated autonomous equation possesses a fixed point will have a background 
solution in the form of a power law [equation (|l2|)l with exponent given by the 
scale invariant condition. 


2.3.2. Linearized solution 

Let us linearize the autonomous equation by setting z(t) = z* + zi(t) + 0[zi(t) 2 ], 
where z* is the fixed point and z\(t) is a small perturbation about z*. F(z(t)) = 0 
reduces to 


F(z(t)) = F(z *) + Fi(zi(t)) + 0[zi(t) 2 ] = Fi(zi(f)) + 0[zi(f) 2 ] = 0. (13) 


So in the linear approximation, we must have 


F 1 (z 1 (t))=0. 


(14) 


But -Fi( 2 i(t)) = 0 is an ODE which is both autonomous and linear, and still re¬ 
tains its character as a nth-order differential equation. Autonomous plus linear is a 
very powerful constraint, and in fact implies a differential equation with constant 
coefficients. That is, there is some polynomial P{- ) such that the function Fi(zi(t)) 
satisfies 


*i(*i(i)) = P 


2l(f). 


(15) 


The standard techniqueB for solving differential equations of this type is to insert a 
trial solution of the form z\{t) = Aexp(Af). Then the differential equation reduces 
to the simple nth-order polynomial 


P{ A) = 0, 


(16) 


which has at most n distinct solutions. Assuming for the time being that all roots 
are distinct, the general solution to the linearized autonomous equation is 


-gift) = Y A i exp(A jt). (17) 

i= 1 


That is 

n 

z(t) = z^ + ^^Ai ex.p(\it) + 0(A 2 ). (18) 

i= 1 

In terms of the original variable y(x) 


y{x) = x p 


n 

z * + 'y ' Ai x + 0(A 2 ) 


(19) 


i—l 

Near the fixed point, we now see the presence of multiple power laws. The pref¬ 
actor x p is governed directly by the scale invariant condition, while the subsidiary 
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exponents arise from solving a polynomial equation based on the linearized approx¬ 
imation. 

We are now in a position to formulate the following theorem: 


Theorem 2: Any nth-order ordinary differential equation that is scale invariant 
and whose associated autonomous equation possesses a fixed point will have a lin¬ 
earized solution which is generically expressible in the form of a power law (equation 
fi9|) in a neighborhood surrounding the fixed point. 


We can also enunciate a general rule of thumb: Because the exponent p ulti¬ 
mately arises from algebraic manipulations involving the physical dimensionalities 
of the various quantities present in the ODE it is likely to be a rational number. In 
contrast the exponents A*, being solutions of a polynomial equation, will generically 
be irrational. 

In the exceptional case where the roots of the polynomial P{ A) are not distinct, 
then there are additional technical complications. Let m of the roots be distinct, 
and let g, with i £ [1, m] be the multiplicity of the zth distinct root. Then we have 
J2T 9i = n and 


m gi — 1 

2i W = EE Aij f 1 exp(A,f). 

i=1 3 =1 


That is 


m gi -1 

z(t) = z* + EE Aij f 1 exp(Aif) + 0(A 2 ). 

*= 1 3 =1 

In terms of the original variable y(x) 


( 20 ) 


( 21 ) 


y{x) = x p 


m gi — 1 

Z. + EE A V ( lna; ) i_1 4 + 0(A 2 ) 

*= 1 3=1 


( 22 ) 


Thus logarithmic deviations from power law behaviour, while not generic, are cer¬ 
tainly relatively easy to get - despite the lore based in statistical field theory and 
second-order phase transitions a logarithm is not necessarily the result of one-loop 
physics. Logarithms can arise from causes as mundane as a repeated root in a 
linearized autonomous ODE. 


2.3.3. Frobenius-like solution 

Going beyond the linearized approximation, in the vicinity of the fixed point, one 
can write the expansion more systematically as z(t) = z* + ezi(t) + e 2 Z 2 (t) + 0[e 3 ], 
where e <C 1, in which case F(z(t)) = 0 reduces to 

F(z(t)) = F(z *) + eFi((zi(t)) + e 2 [Fi(^(i)) + F 2 (zi(i), zi(t))] +0(e 3 )- (23) 
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Here F 2 (■, •) represents a second-order functional derivative of F(-) around the fixed 
point, which will be some messy quadratic differential operator with constant coef¬ 
ficients. This implies both 


F 1 (z 1 (t)) = 0, (24) 

and 

Fi(z 2 (t)) = -F 2 (zi(t),z 1 (t)). (25) 

That is 

- 22 (t) = -Fi -1 [F 2 (zi(t),zi(t))] . (26) 

But we already know that zi = X^"=i A * exp(Aji), so z 2 consists of pieces of the 
form exp([Aj + A j]t). Indeed 

n n 

z 2 (t) = ^ E B ^ Ai ’ A J') eX P(t A i + A #)- ( 27 ) 

»= 1 1=1 


with B(Ai,Aj ) being some complicated function of the A,; that is deducible from 
the precise form of F\ and F 2 . 

Collecting terms 


z(t) — -f 'y ' Ai exp(A^t) + EE B(Ai, Aj ) exp([Ai -I- A ,j]t) + 0(e 3 ). (28) 

2 =1 2=1 j=l 

In terms of the original variable y(x) 


y(x) = x p ■ 


n n n 

z,+EAi Ai +EE B ( A ^ A o) xXi+Xj + °( e3 ) 

2=1 2=1 j= 1 


(29) 


The general pattern is now clear. Let us take A = (Ai, A 2 ,... A„) and let m = 
(mi, m 2 ,.. • m n ) € IV™. Then formally 

z(f) = ^ AL(m) exp ^A ■ fht'j . (30) 

mGW n 


This series may be asymptotic, in general we make no claims about convergence, 
though in specific cases it often does converge in an open region surrounding the 
fixed point. The coefficients A(in) are certainly not independent but are inter¬ 
related via the underlying differential equation. This is basically Liapunov’s ex¬ 
pansion theorem, as discussed for instance by Lefschetzcl. Readers interested in 
additional mathematical rigour are directed to that reference. 

Near the fixed point our equation becomes, in terms of the original variables 

y{x), 

y(x) = x p ■ ^ A(rh) m 

m£iV n 


( 31 ) 
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This type of series is Frobenius-like in the sense that the prefactor x p will be 
governed by the background solution, while the exponents of the power series will 
be given by the linearized problem. Hence, we can formulate the following theorem: 


Theorem 3: Any nth-order ordinary differential equation that is scale invariant 
and whose associated autonomous equation possesses a fixed point will generically 
have a formal solution in the form of a Frobenius-like power series (equation ©) 
in a neighborhood close to the fixed point. 

As per the previous discussion, this form holds only if the roots of the linearized 
problem are distinct. For coincident roots one should expect the Frobenius-like se¬ 
ries above to be modified by various logarithmic terms. Indeed based on Liapunov’s 
expansion theorems the coefficients A(m) become polynomial in t, and so polyno¬ 
mial in In z [cf. equation (^2|) above]. The resulting series is too clumsy to be of 
immediate use and fortunately none of the specific examples we encounter exhibit 
this particular behaviour. 


2.4. Systems of nth-order scale-invariant differential equations 

2.4.1. Background solution 

All of the previous analysis can be generalized to systems of I coupled nth-order 

ordinary differential equations. There are three possible routes: 

(1) Develop a suitable “vector” and “matrix” notation. 

(2) Use the fact that a system of / coupled nth-order differential equations can 
always be reduced to a system of N = I x n coupled first-order differential 
equationsEJ, but you would still need to develop a suitable “vector” and “matrix” 
notation. 

(3) Use the fact that a system of / coupled nth-order differential equations can 
generically be “decoupled” by being reduced to a single DE of order N = I x n 
in one of the dependent variables, though there seems to be no constructive 
algorithm for doing so.tl 


On balance, we have found it most useful to directly develop a suitable “vector” 
and “matrix” notation. Let us begin by considering the equation 


[F] = 0 (32) 

where [A] is the column vector of ODEs 


( F1 

(x,yi(x),y[(x),. 

.,y 2 (x),y' 2 (x),. 

■ ,yi(x),.. 

d n yi(x) 

* ’ dx n 

f 2 

yx,yi{x),yi(x),. 

. ,y 2 {x),y 2 (x),. 

-,yi{x),-- 

d n yi(x) 
’ dx n 


(x,yi(x),y'i{x),. 

■ ,V 2 (x), y' 2 (x),. 

-,yi(x),-- 

d n yi(x) 
* ’ dx n 


( 33 ) 
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When we say that this system is scale invariant, we mean that as x —> ax and 


where 


[r(*)]^[A][r(*)], 


[A] 


~ a Pl 0 0 ... 0 " 
0 a P2 0 ... 0 

0 0 0 ... a? 1 


( 34 ) 


(35) 


and [Y(a;)] is a column vector of dependent variables, the equation [F] = 0 re¬ 
mains invariant. Note that this means that we have chosen our dependent variables 
Ui(x) ... yi(x) in some convenient manner to make their scale transformation prop¬ 
erties both simple and independent of each other - this will almost always occur 
automatically in systems of ODEs derived from an underlying physical problem. 

Previously, we also showed that any nth-order differential equation in a single 
dependent variable that is scale invariant can be put into its associated autonomous 
form by applying the substitutions described in equations (^) and <©■ Therefore, 
it carries over that we can make the same kind of substitution in each equation of 
the system. After doing so, we can define a new column vector [F] in the following 
manner: 



^ Fi ^2i(f),ii(t),.. 

.,z 2 (t),z 2 (t),. 

.,zi(t),.. 

d n zi(t) ^ 
‘ > di n , 

\ 

[F] = 

F, ^i(i),z'i(t),.. 

-,Z2(t),Z 2 (t), ■ 

-,2/(f),.. 

d n *i(tr 

* ’ dt n j 




.,z 2 (t),z 2 (t),. 

■,Zl(t),.. 

d n zi(tF 

’ ’ d t n j 

J 


(36) 


where for each i £ [1,7] we have Di(x) = x Pi Wi(x) = e tpi Zi(t) and furthermore 
x = e t . From the previous analysis we know that each equation in this new vector 
equation [F] = 0 is now autonomous, hence, the whole system is autonomous. 
Let’s now assume that this autonomous system has a fixed point [Z*] such that all 
derivatives of z(t) with respect to t in [F] vanish at the fixed point, where we mean 


fzt\ 



\zV 


(37) 


[Z(t)\ = 

( zi(t)\ 
z 2 (t) 

= [Z*] = 

(zl\ 

Z 2 


\Zl{t)j 


\z*J 


( 38 ) 


That is 
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is one specific solution to the system of ODEs. 

Substituting back we obtain the following equations for [Y(a;)]: 


/ yi(x)\ 


■e 

o 

o 

o 


( Z *A 


/ Z* X P1 \ 

V2(x) 


0 x P2 0 ... 0 


Z 2 


Z% X P2 

Kyi{x)j 


1 

.. a 

• • o 

• • o 

• • o 

_ 1 




\ z *J X PI ) 


Therefore, we can formulate the following theorem: 


Theorem 4: Any system of I nth-order differential equations that is scale invari¬ 
ant, and whose associated system of autonomous equations possesses a fixed point, 
will have solutions in the form of a collection of power laws (equation |39| ), with 
exponents given by the scale invariant condition. 


2.4.2. Linearized solution 

If we linearize the nth order system, then for each F % we will obtain a linear differ¬ 
ential equation that has the capability of mixing the various Z{. Specifically let 

Ziit) = z* + ez lt i(t) + 0(e 2 ), (40) 

then 

n 

(*)) = Fi(^) + + 0(e 2 )- (41) 

3 = 1 

Here the { F\ }are linear differential operators that lead to a system of linear and 
autonomous differential equations 

n 

£{A}y(*i ,j(t))=0. (42) 

3=1 

This implies that each of the {Fi}ij is a constant coefficient polynomial in d/dt 

mu = p i3 (^) • ( 43 ) 

Inserting a trial solution of the form Zij(t) = Aj exp(At) now yields the vector 
equation 

n 

X)A,=0. (44) 

3=1 

There is a non-trivial solution if and only if the I x / matrix Pij(X) is singular 
(det/Py/A)} = 0), in which case the Aj correspond to the singular eigenvector. 
Since each element of the matrix Pij{ A) is itself a polynomial of order n, (where n 
is the order of the individual differential equations in the system), the determinant is 
a polynomial of order N = Ixn having up to N separate roots. Let us assume these 
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roots are distinct and call them Aq, with a G (l,iV) and denote the corresponding 
eigenvector by (A a )i. As long as the roots are distinct these eigenvectors, multiplied 
by the associated exponentials, will span the solution space and we can write the 
general solution to our autonomous system as 

N 

= ^2(A a )i exp(A a t). (45) 

a=l 


In terms of our original equation, 


N 


Ui(x) = x Pi < z* + ^2(A a )i x Xa + 0(e 2 ) > . 


(46) 


a=1 


This now allows us to formulate the following theorem: 


Theorem 5: Any system of I nth-order differential equations that is scale in¬ 
variant and whose associated autonomous equation possesses a fixed point will 
generically have solutions in the form of a collection of power laws (equation [46]) 
in a neighborhood close to the fixed point. The prefactors \x v ] will be governed by 
the background solution, while the X a are given by the linearized problem. 


If the roots are not distinct then a much messier equation involving logarithms 
in x can be constructed, but is beyond the scope of the applications we have in 
mind. 


2.4.3. Frobenius-like solution 


We can again go beyond the linearized approximation in the vicinity of the fixed 
point and write Zi(t) = z* + ez\^(t) + e 2 Z 2 ,i(t) + 0[e 3 ]. Then, 


Fi(zj(t)) = Fi (z *) +eJ2{F 1 } ij {z ltJ ) 

3 =1 


+e 


2 


N 


N N 


^2{Fi}ij(z 2 ,j) + Zi,k(t)) 

j =i j=i fc=i 


+ 0(e 3 ), 


(47) 


where, similarly to the preceding discussion for a single dependent variable, the 
quantity {-FbjijfcO, 0 represents the second-order functional derivative of {F}j(-) 
around the fixed point [Z*]. Using the same reasoning as before, we will obtain a 
solution of the form 


N N 

Z2,i(t) = Y Y{ B ( Aa ’ ex p([^« + (48) 

a=l (3=1 

Collecting terms and substituting back to our original solutions of yi(x) we obtain 
N N N 

Vi(x) = X Pi Z* + + Y Y( B ( A ^ A ^ xXa+XP + °( £3 ) | > ( 49 ) 

a=l ot =1 (3=1 
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The general pattern is again clear 

Vi{x)=x Pi ■ Y, (50) 

m£]N N 

where as before, A = (Ai, A 2 ,... A n) and m = (mi, m 2 , ■ ■ • m^v) £ 1N N . We can now 
formulate the following theorem: 

Theorem 6: Any system of / nth-order differential equations that is scale in¬ 
variant and whose associated autonomous equation possesses fixed points will 
generically have a solution in the form of a Frobenius-like power series (equation 
[To| j in a neighborhood close to the fixed point. The prefactors x Pi will be governed 
by the background solution, while the exponents of the power series will be given 
by the linearized problem. 

As previously, there is an explicit requirement that the exponents arising from the 
linearized problem be distinct, otherwise (as per Liapunov’s expansiormtheorem) 
messy expressions polynomial in logo; will be encountered. (See Lefsclietzcl for tech¬ 
nical details.) 


2.5. Discussion 

Before we turn to several specific astrophysical examples to show how this formalism 
works in detail, let us summarize key points: 


Scale invariance does not automatically imply power-law behaviour; the ex¬ 
istence of one or more fixed points is an additional requirement. An explicit 
example of a scale-invariant ODE without any power-law solution is presented 


in Appendix A . 

If we perturb slightly away from the fixed points, there will no longer be an 
exact power-law behaviour, but we have something similar: Generically there 
will be a Frobenius-like expansion of the form 


OO 


8=0 


y(x) = x e ■ y ai x 


(51) 


where the exponent p depends on the fixed point and the exponent A depends 
on the linearized ODE around the fixed point. The coefficients at are of course 
constrained by the original exact ODE, and since it is nth-order we expect n 
of these coefficients to be “free” (that is, to depend on initial conditions) with 
the other coefficients in principle being determinable from the first n and the 
exact ODE. 

• We feel that this formalism, attributing as it does power-law behaviour to rather 
general features of differential equations, goes a long way towards explaining 
the apparent ubiquitous occurrence of power-laws in nature. 
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3. Application: Isothermal Newtonian stars 

In the following section, we apply the general formalism developed above to a 
static isothermal Newtonian star. The relevant system of ODEs is scale invariant, 
and contains two fixed points, one corresponding to the center of the star, and the 
other corresponding to the point at infinity (where the star develops an infinitely 
thin halo). We formally solve the differential system associated with it and obtain 
power-series solutions for the total pressure, radius, compactness and mass of the 
star in terms of a Frobenius-like series about either fixed point. For each of these 
physical quantities, the radius of convergence of these series (one for each fixed 
point) overlap, so they constitute a complete power-series solution. These results 
subsume and extend previous results found by Chandrasekhara, and the classic 
results of Lane and Emden. In particular, these series are then used to discuss the 
behaviour of such stars near the onset of collapse. 


3.1. Background solution 

The basic equations of (static) non-relativistic stellar structured are (G = 1) 

d p m(r) 

t = ~p'< r > — ~; 

dr 


i(r ) = f 4-7T r' 2 p(r') dr / . 

Jo 


(52) 

(53) 


Let us now adopt the following particularly simple equation of state, which is ap¬ 
propriate, for instance, for an isothermal star: 


P(r) = c 2 s p. 

(54) 

Then 


d p p(r) m{r) 

dr c 2 r 2 

(55) 

dm(r) 47rp(r)r 2 

dr c 2 

(56) 

where c s stands for the speed of sound. The relevant boundary conditions are 
m(0) = 0 and p(0) = po- This system of differential equations is scale invariant 
under r —> ar, p —> a~ 2 p and m —> a +1 m. This means that from our previous 
analysis we have (j/i, 2 / 2 ) = (p,m) and (^ 1 ,^ 2 ) = (—2,1). As stated previously, we 
can make suitable transformations to convert this scale invariant system into an 
autonomous system. For this purpose, let’s substitute 

C(r) = 27rp(r)r 2 /Cg, 

(57) 

and 


X(r) = m(r)/ (2 c 2 s r) . 

(58) 
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Then, the system becomes 


r ^ = 2C(r) [1 - X (r)]; 


dr 

„ d.x( r ) 

dr 


= C (r) - xir). 


(59) 


(60) 


As expected this system is now equidimensional-in-r. Define a new variable t = ln(r) 
so d/di= rd/dr, then 

dC (t) 


d t 

dxft) 

d£ 


= 2 C (t) [i - x(t)]; 


= CO) - x0)- 


(61) 


(62) 


As we can clearly see this system is now autonomous. There is a fixed point at 
X = 1 and £ = 1. For any initial condition this point is an attractor as t —» oo. The 
center of the star is at r = 0 (f = —oo). For a regular star with finite po we have 
( = 0 and x = 0 as t — oo; which is the second fixed point. So solving our special 
dimensionless differential equation is equivalent to looking for the unique integral 
curve that emerges from the fixed point at (0, 0) and terminates on the fixed point 
at (1,1). Applying Theorem ^ we can easily check that a singular solution to the 
system is of the form 


p(r) = C\ r 


-2 


(63) 


and 


m(r) = C 2 r 1 , (64) 

where C\ and C '2 are constants. Solving for the constants we see that the (1,1) 
fixed point corresponds to the power-law solution 


and 


i(r) = 2 1 


p{r) = 


2nr 2 ’ 


(65) 

( 66 ) 


while the (0, 0) fixed point corresponds to the more prosaic p(r) = 0 = m(r). 


3.2. Power series in terms of the radius 
3.2.1. Linearization about the fixed points 

In order to obtain the linearized solution it is necessary to linearize about its critical 
points. Let us first make the substitution x = Xi + 1 an d C = Ci + 1- Then, 

dCi(t) 


d t 

d\i(t) 

d£ 


= —2xi (t); 


= Ci(t) -xi(*)- 


(67) 

( 68 ) 
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Close enough to the (1,1) stagnation point, solutions are of the form 



0 —2 
1 -1 


(t ~ U ) 


\x(u)J 


(69) 


The eigenvalues of this matrix are 5 = (— l±iy/7) /2, and the approximate linearized 
solutions approach the critical point exponentially in t — as exp(<5 1) = r s . 

There is a second (unstable) stagnation point at (0,0) — a repeller. If we lin¬ 
earize around ( 0 , 0 ) we get 

^ = +2Ci(i); (70) 

^^=Ci(t)-Xi(t). (71) 

The corresponding eigenvalues are ( 71 , 72 ) = (2,-1), corresponding to regular r 2 
behaviour for f and x a t the center of the star, plus an otherwise possible but 
irregular and unphysical behaviour for \ that goes as 1/r. This 1/r behaviour 
would correspond to a finite m( 0 ) if such were possible. 


3.2.2. Power series about (0,0) in terms of the radius 

Based on Theorem and our general comments regarding Frobenius-like series it 
is reasonable to assume the following expansion to solve the system 

OO OO 

xw = EE i - zra71+n72 - ( ?2 ) 

m= 0n—0 

Here we have introduced the dimensionless variable 

z = \J 4tt p 0 /c 2 s r = (y/^po/cl) r, (73) 


with the particular numerical coefficients being chosen for future convenience. In 
this case, since the 72 = — 1 index is unphysical, it will be excluded from the 
analysis. Hence, we restrict attention to an expansion of the following form, 

OO OO 

x( z ) = J2 b - zUJ1 -* E- 2n - ( 74 ) 

n —0 n—0 

In order to solve for the correct coefficients b n we need to plug the previous equa¬ 
tion into our differential equations. First, by uncoupling the system (eliminating f 
to obtain a single second-order differential equation for \) we obtain 

d 2 x(z) __ 2 x (z) 2 X (z) 2 2 X {z) dx(g) _ 

i 0 2 ' 2 ' I V ^) 

az z z z z z z d z 

The relevant boundary conditions may be inferred from the fact that the pressure 
(and so also the density) is finite at the origin: 

x(0) = x'(0) = 0; x"(0) = i; 


X' = dx/dz. 


(76) 
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The fact that only integer exponents arise in our series solution for x( z ) is ultimately 
due to the fact that the second-order ODE above, though nonlinear, has what 
amounts to a regular fixed point at z = 0. So standard results in terms of ordinary 
Frobenius series suffice to develop an appropriate ansatz. 

Now multiplying by z 1 , we can write 

2 , 1 


dz 


z~T X - z ~rX - 2x = - 2 % z—x, 


dz 


or in operator form, 


dz 


z— - z- - 2 


dz 


X = 


2 + 4 

dz 


dz' 


X 


Now we plug in our assumption of the expansion in equation 
recursion relation of the form 

n 

[( 2 n ) 2 - ( 2 n) - 2 ] b n = - [2 + ( 2 n)] ^ &i b n - x . 


i =0 


Because 60 = 0 we can rewrite this as 

(n + l) 


bn — 


(n — l)( 2 n + 1 ) 


^ ' bi b n -i- 


(77) 

(78) 

2 ) and we get a 

(79) 

(80) 


i -1 


This recursion relation, together with the initial condition b\ = 1/6, derived from 
the fact that x"(0) = 1/3, completely specifies the power series about the (0,0) 
critical point. It is straightforward to compute these coefficients numerically, though 
there does not seem to be any simple explicit closed form for the general term. 

For £ we use the relation 


Z dz X = C - X- 


So if 


then 


CM = £ 


a n z 


n —0 


= (2 n + 1) bn 


(81) 


(82) 


(83) 


Keeping the first 250 terms in this series, and comparing with an explicit 4th order 
Runge-Kutta integration, gives very high accuracy [O(10 -6 )] out to z « 3. 

In terms of the physical variables m{r) and p(r ) we have 


oo / . 2 \ n 00 

ro(r) = 2 c 2 s rx(r) = 2 c 2 s r ^2b n ( ) = 2 c 2 s r ^ b. 


n= 1 


f 47rp(p 


n =1 


(84) 
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and 


p(r) = 


CO 


2irr 2 


2ixr 2 

71=1 


/ 4np 0 r 
\ c 2 


2irr 2 ' 

71=1 


/ 4np 0 r 

\ cf 


(85) 


These are Taylor series for the solution to the differential equations that describe 
an isothermal star, with the coefficients being specified by the recursion relations 
determined above. By bounding the coefficients using the recursion relation, one can 
show that the series definitely converges for \z\ <3 and that it definitely diverges 
for \z\ > v/12. Numerical methods seem to indicate a radius of convergence of 
R ss 3.273687274. The a n and b n coefficients do not have a simple explicit analytic 
form, but computing them numerically from the recursion relation is trivial. 

Note that in more classical language this analysis is equivalent to finding a power 


series solution of the isothermal Lane-Emden equation. See Appendix C for details. 


3.2.3. Power series about (1,1) in terms of the radius 

Turning to the other critical point, from equation ( 0 ) we see that x(°°) = 1- Now 
making a shift of variables appropriate to the (1,1) fixed point, \ = 1 + X'i> where 
Xi is small, we obtain 


( d d 

\ Z dz Z dz 




( 86 ) 


with \i(oo) = 0. Let’s assume an expansion for xi of the form given in equation 
(0) with the substitution ( 71 , 72 ) —► (<5i, £> 2 ). The real parts of < 5 y 2 must be negative 
so that the expansion has the right limit at infinity, where xi —> 0. By construction 
we know that boo = 0. In order for x to be real, the matrix of coefficients b mn must 
be Hermitian. Now let’s insert this ansatz into our shifted differential equation 
[equation (|| 6 [)]. For convenience define A = mSi 4 - nf> 2 , since then 

z ^_ z m5 1+ n5 2 = Z _± Z A = A ^ ( 87 ) 

d z d z 

and so equation (^) becomes 

77i n 

(A" + A + 2) b mn = — (A + 2) bij b m _i- <n _j. ( 88 ) 

2=0 j= 0 


This implies the recursion relation 


b 


mn 


(A+ 2) 

(A 2 + A + 2) 


'Y, 'Y, kj b m -i- n -j. 

2=0 j =0 


(89) 


Since boo =0, the right hand side does not involve b mn itself. Thus b mn can be cal¬ 
culated once the lower-order bij are known. In particular once 601 = is specified, 
all other coefficients can be calculated recursively. Unfortunately the only known 
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way of calculating &oi = ^io is by performing a numerical integration ( e.g ., Runge- 
Kutta) out to large radius and then numerically fitting the data. Summing up, we 
obtain 

OO OO 

X (z) = 1 + £ £ b mn . (90) 

m=0 n =0 

A similar power series will exist for £(z). Indeed suppose 

OO OO 

C(z) = 1 z m5 ' +nS *. (91) 

m= 0n—0 

Then 


Q"mn — T 1) bmn — {jTld\ T T 1 S 2 T 1) b mn 


(92) 


In terms of the physical variables m{r) and p{r) 

( OO OO 

m(r) = 2 c 2 s rx(r) = 2 c 2 s r < 1 + ^ b mn 

and 


m —0 n=0 


/ 47T p 0 r* 

\ 4 


2 \ (m<5i+n<52) 'j 


p{r) 


C (z)ct 

2irr 2 



(93) 


(94) 


These are Frobenius-like power series with fractional complex exponents for the 
isothermal star [recall <5^2 = (—1 ±*\/7)/2], where the coefficients a mn and b mn are 
specified by the recursion relations determined above. These expansions are con¬ 
siderably more general than the discussion presented in Chandrasekhar 13, and the 
related results for the relativistic case implicit in Harrison et aZ.i. Those results are 
tantamount to just keeping the first sub-asymptotic term, which is also contained 
in our solution, namely 


and 


i(r) 


2c 2 r ■ 


1 + 23? 


no 


/ 47rp 0 r 


p{r) 


2'kv 2 


1 + 23? 


aio 


( Attpqt 



(95) 


(96) 


Using analytic techniques, the series can be shown to converge for |z| > I|b|5oi| 2 . 
This is a weak bound but at least establishes a non-zero radius of convergence. Nu¬ 
merically fitting | 6 qi | by matching to a 4th-order Runge-Kutta integration implies 
guaranteed convergence for |^| > 5.246853637. Direct numerical evidence suggests 
that the power series is reasonably well behaved down to z ~ 3 with the first three 
terms providing 0.5% accuracy. By numerically applying the ratio test, we find 
that this power series seems to converge for \z\ > 0.165827. Indeed if we keep the 
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first 100 “levels” {n + m < 100) then the power series reproduces a Runge-Kutta 
integration to within 6 significant figures for 2 > 0.21. 

The combination of both power series, a Taylor series around r = 0 and a 
fractional power series around r = oo, provides a solution everywhere to arbitrarily 
good accuracy, since the two expansions have overlapping radii of convergence. 


3.3. Power series in terms of the pressure 


3.3.1. Linearization about the fixed points 

Instead of expanding x and C in terms of radius r, we now want to develop an 
expansion in terms of pressure p. This will be useful if we want to generalize beyond 
a simple linear equation of state. 

We can rewrite the equation of hydrostatic equilibrium in the following manner: 

dp 2 p(r)x(r) 

dr r 

Now, by applying the chain rule to d£/dp, differentiating equation ( |59| ) with respect 
to r, and plugging in our new expression for the equation of hydrostatic equilibrium, 
we obtain 


(97) 


dC , C 

P t = C-• 

dp X 


(98) 


The same can be done to equation 

dx 


P 


= - 1 - 


to obtain 

C' 


(99) 


_ 1 

dp 2 v X, 

This new system of differential equations in terms of the pressure still has a stagna¬ 
tion point at (C> x) = (1)1)- This is the same stagnation point that we encountered 
earlier. This stagnation point is again an attractor, which means that for any initial 
condition as p —> 0 almost all solutions of the system approach this solution. Note 
that the second stagnation point (the repeller) at (0,0) is still there but somewhat 
disguised, since the system in terms of p is highly singular there. We will not need 
the second stagnation point in the following discussion. Linearizing about the (1,1) 
stagnation point we obtain 


dCi 

p -j— = xi; 

d p 

d Xi 1 , , x 

P W = 2 (X '- <1) - 


( 100 ) 


( 101 ) 


Close to the (1,1) stagnation point, the solutions to the system are of the form 


( C 'f A - 

VXi(p). 


exp 


In ip/Pi) 


( Ci (pi. 
\XiiPi 


( 102 ) 


0 1 

— 1/2 1/2 

The eigenvalues of this matrix are A = (1 ± iy/l)/A. Note that this is simply —S/2, 
where S are the exponents as encountered when working with the variable r. 
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3.3.2. Power series about (1,1) in terms of the pressure 


It is essential to decouple the system of differential equations in order to obtain a 
power series expression for the mass and the radius in terms of the pressure That 
is, we want to replace the first-order system by a single second-order equation. In 
order to decouple equations ( fisl ) and (^9j), we solve for £ in terms of x from equation 


C = X 1 - 2 p 


,dx 

dp 


(103) 


Now, if we insert this expression back into the first differential equation, we obtain 


- 1 



XP 


,dx 

dp 


P 


dx 

dp 


x p-t- p~r x 

dp dp 


(104) 


We can see from the analysis in the previous section that \{r = oo) = \{p = 0) = 1. 
It is convenient to write x = 1 + X because we now have \{p = 0) = 0. After this 
shift in variables, our decoupled differential equation becomes 

d d 


d d cl 

1 _ P~T + 2 PT~ Pa~ 

dp dp dp 


X = 


d,X 

P dp P dp P dp 


x 2 . 


(105) 


Now, in line with Theorem g and our general discussion regarding Frobenius-like 
series, we assume an expansion of the form 

oo oo / \ mAi+nA2 

x{p) = EE dmn ) ■ (106) 

m=0 n=0 \P0 / 

Previously we used the dimensionless variable z = (y/4frpo/c^)r to create a power 
series expansion about the stagnation point (1,1) for \ an d (■ This time we find 
it more convenient to expand in powers of p/po- It is interesting to note that our 
eigenvalues, Ay 2 , make physical sense because we need positive real parts in order to 
satisfy the boundary condition \(0) = 0. For \ we have doo = 0 and the Hermiticity 
condition discussed previously still holds. Again, let us for simplicity define A = 
toAi + n \2 ■ We then have 


P- 


P_ 

dp \p 0 
And therefore equation 

dmn = 


raAi+nA2 


dp 


= A 


= P 

05|) implies 

A (1 — A) . . 

(1 — A + 2A 2 ) 2^ 2^ *3 

V ’ 2=0 j =0 


(107) 


(108) 


Note that by construction doo = 0, so once doi = d* 0 is specified all other coefficients 
can be calculated recursively. Now this recursion relation holds for the coefficients 
in x only. For \ itself we have 

\ mAi+nA2 

p 


x(p/po) = 1 + E E dmn ( — 

m=0n=0 X 


(109) 
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A similar power series will exist for C( p/po ) = Cip/Po) — 1- Let us assume the 
following expansion 


OO OO 


Cip/Po) = EE' 


P 0 


m =0 n —0 

Solving for from equation (^8|) we obtain 

>■ dy 2 

C = X-P-r -, 

dp 

which, upon shifting to the fixed point, is equal to 


( 110 ) 


( 111 ) 


C = X ~ P ^ (2X + X 2 ) ■ 


Comparing exponents we obtain 


( 112 ) 


m n 

Cmn — ^mn A ( 2(i mn H - EE dj,j d r n—i: 

i=0 j=0 


n-j 


Using the recursion relation for the 

A (2A 2 -A + l) ^ N 

Cmn dmn A I "L A~(l A) ^mn 


2A 

W~A 


(113) 


(114) 


Note that by construction coo = 0. Now this relation holds for the coefficients in £ 
only. For £ itself we have 

oo oo / \ mAi+nA2 

CGVpo) = i + EE c ™ £ ■ (ns) 

m=0n=0 \P / 

We could of course unwrap the definitions of £ and x to obtain Frobenius-like power 
series for the physical observables r(p) and m(p). 


3.4. Mass, pressure, and radius oscillations of a star on the verge 
of collapse (po —► oo) 

Let us now consider a more realistic equation of state: 


Pip) 


( Pc = Pc/c 2 s ■■ P<Pc 

\p/c 2 s : p>p c 


(116) 


where p c and p c stand for the pressure and density at the surface of the core. Most 
of the analysis carried until now still holds, at least for the central region of the star. 
But we now chose to replace the previous equation of state [p(p) = p/c 2 , which was 
linear out to arbitrarily low density] with a new equation that is somewhat more 
realistic. We have a dense core described by a linear equation of state surrounded 
by an envelope of constant density in which the pressure drops rapidly to zero. 
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Since in the previous section we have established the existence of a power series 
for C(p/po), we can formally solve equation © for the radius of the core itself by 
setting p —> p c (at this stage making no comment about the envelope). Hence, 


T core {Pc: PO ) 



\/C(Pc/P0) 


c 


2 

s 


\Z2np c \ 


OO OO 

i + EE C-mn 
m= 0n—0 



(117) 


This last equation clearly shows the oscillatory behaviour of the core radius as a 
function of central pressure, due to the imaginary exponents in the power series. 

If we now look at equation (|58| ) and solve for the mass of the core in terms of 
the compactness and radius we obtain 


m coie (pc,po) = 2 c 2 s T core (Pc,Po) X'(Pc/Po) 

2 ci 


\Z2np c ^ 


OO OO 


1 + EE< 

m= 0 n =0 


OO OO 


1 + E E 


m =0 n —0 


(118) 

A') 


As we can see from this equation the mass of the core will also oscillate as a function 
of central pressure due to the imaginary part of the exponent of the power series. 
Indeed there is a related power series such that 


TOcore(PcPo) = 


2cf 


\/27 rpc 


i + EE’ 

k m=0 n —0 


(119) 


with the Hermitian coefficients m mn being in principle calculable in terms of the 
Cmn and (Iran • 

The above expression gives us the mass of the core of the star, not of the star 
itself. In order to calculate the total mass of the star we need to first calculate the 
mass of the envelope of the star: 


WenvelopeM = / 4:TTr' 2 p c dr . 

d r mrfi 


( 120 ) 


Hence, 


mtotai(r-) = -TT Pc (r 3 - r 3 ore ) + : 


( 121 ) 


We calculate r sur f ace from our equation of state and equation of hydrostatic equi¬ 
librium: 


Penvelope(r) = Pc ~ f Pc dr', (122) 

Acore 

where we have used the fact that p —> p c is constant throughout the envelope. 
Integrating 

PenvelopeM = Pc + (^Pc™ core — 7^ 7r Pc r core^ ~~ ~ ^ A Tj^Pc ( 7 ' core _ ^ ) ‘ 
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We can solve for r sur face by solving for p e nveio P e(r = Surface) = o. The resulting 
cubic equation is 

Penvelope (c = r sur f ace ) = 0 = C\ r gur £ ace + C 2 ^surface O 3 , (124) 


where we have chosen the constants C\, C 2 and C 3 to be 

2 9 

Cl = 


C 2 = 


p c m c 


—^"/+ Core Pc, 


c 3 = p c m cole - -7T p 2 c r] 


, 2 r 3 

c core' 


(125) 

(126) 

(127) 


This cubic equation can be explicitly [if tediously] solved for r sur face as a function 
of p c , m C ore, and r cor e- Suffice it to say that we now have implicitly solved for the 
total mass of the star (M tota i), where we have taken r sur f ace as a parameter (which 
can ultimately be expressed, via m core , and r core , as an oscillating function of p c 
and p 0 ). That is 


d7total(Po Po) "Al (^"surface (j?c: PO ) I'core (pc. ? Vi) ) )T H^c.ore (Pc ■ Pi) ) ■ ( 1 28) 

Ultimately, since the individual terms above exhibit oscillatory behaviour, there 
will be some Frobenius-like power series such that 


AAotal(PcPo) = A/total (j>cP0 00) 


1 + EE m - 

„ m =0 n —0 



(129) 


From this expression for the total mass of the star we can clearly see the oscillatory 
behaviour mentioned earlier. 

Note that this behaviour is much more extensive than the scaling results dis¬ 
cussed in Chandrasekhar 0 [see also Harrison et ala \ . Their result is equivalent to 
keeping only the first nontrivial term in the power series. That is: 

Af total (p c ,p 0 ) « -Mtotai(PoPo -> oo) <1+23? 


-Mm [ — 
Po 


Si 


(130) 


We now see that the scaling behaviour they discussed is only the first term in an 
infinite series. 

It may be objected that our current analysis is only valid for the special case of 
a linear equation of state with a constant density envelope. However, from the way 
we have set up the general theory it is clear that this restriction is not particularly 
serious. First, any reasonable equation of state will saturate to a linear at high 
enough densities, so our results can be applied without modification to the core of 
any reasonable stellar configuration as the central pressure po —> oo. Furthermore, 
suppose one has some generic equation of state. Then the exact system will not be 
scale invariant. Nevertheless there will be a special critical configuration in which 
the central pressure is infinite, and since by hypothesis the equation of state becomes 
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linear at sufficiently high pressure, the critical solution will exhibit approximate 
scale invariance in the high-pressure region — and this will correspond to a fixed 
point in the autonomous system. Linearization around the fixed point will generate 
some critical exponents which depend only on the high-pressure limit of 

the equation of state. But deviations from the scale invariant critical configuration 
in the core of the star will now feed into the exact scale non-invariant equations, 
still leading to the total mass being given in terms of some Frobenius-like series. 
Ultimately this Frobenius-like power series behaviour is generic to many systems 
of differential equations and intrinsically does not have anything to do with gravity 
per se. 


4. Application: Isothermal relativistic stars 
4.1. Background behaviour 

In order to treat relativistic stars, we need to consider the relativistic equation of 
hydrostatic equilibrium (TOY equation) 


dp _ / p\ (to + 47rr 3 p/c 2 ) 

dr v c 2 / r 2 (1 — 2 m/(rc 2 )) 


(131) 


We again have a system of two first-order differential equations, namely equations 
(131) and the mass equation (|56|), which must be supplemented by an algebraic 
equation of state. Note that in the TOV equation c is the speed of light, and that 
as c —> oo we formally recover Newtonian physics. We will consider an equation of 
state identical to that used in the previous chapter, i.e. p(r) = c 2 s p(r), where c s 
stands for the speed of sound in the fluid. Making the appropriate substitutions we 
obtain the closed system 


S 



< 

dm 4irr 2 p 

dr cl 
* 


(to + 47rr 3 p/c 2 ) 
r 2 (1 — 2 m/(rc 2 )) ’ 


(132) 


Note that (as in the Newtonian case) this system S is invariant under the substitu¬ 
tions r —> ar, p —* a~ 2 p and m —> a m. Revisiting our theoretical analysis we also 
discover that, ( 3 / 1 , 2 / 2 ) = (p, to) and (pi,P 2 ) = (—2,1), again as in the Newtonian 
case. 

It is possible to systematically transform the system S into an equidimensional- 
in-r system Si, and subsequently into an autonomous system S 2 . Let us first 
substitute ( = 2np(r)r 2 /c^ and x = m/(2rc 2 ). Let us also define t = ln(r) so 
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d/d t = rd/dr. Then, the system becomes 


Si 


dC 

d t 


2C^2C(l + /3 2 ) 


(x + P 2 C) 

1 - 4/3 2 x ’ 


dx 

dt 


= C-X) 


(133) 


where /? = c s /c. As we can clearly see S/ is now autonomous, and therefore, 
translation invariant. There are two fixed points, namely (a, a) and (0,0), where 
a = 1/(1 + 6/3 2 + (3 4 ). Note that /5 —> 0 corresponds to the Newtonian limit. 

As we found previously, the first fixed point is an attractor as t —> oo, which 
corresponds to 1/r 2 behaviour for the pressure. The second fixed point at (0, 0) is 
a repeller, as we also found in the Newtonian case. This fixed point corresponds to 
regular r 2 behaviour for £ and \ at the center of the star, plus an irregular and 
unphysical behaviour for \ that goes as 1/r. The center of the star is at r = 0 
(t = —oo). For a regular star with finite po we have £ = 0 and x = 0 as t —> —oo. So 
solving Si is equivalent to looking for the unique integral curve that emerges from 
(0,0) and terminates at (a, a). Since system S is scale invariant and its associated 
autonomous system Si has fixed points, we can apply Theorem [| and easily check 
that the background solution to S is 


P(r) = Ci r 2 , 

and 

m(r) = C 2 r 1 , 


(134) 


(135) 


where C\ and Ci are constants. Inserting back into S and solving for the constants 
we see that the (a, a) fixed point corresponds to the power-law solution 


and 


m(r) = 


2c 2 


1 + 6/3 2 + p 


■2 , «4 ’ 


p(r) = 


1 


2-7rr 2 1 + 6/3 2 + /3 4 


(136) 


(137) 


Note that if we take the limit c —> oo, we recover the exact solutions obtained in the 
previous chapter for the Newtonian case. Also note that the background solution 
for the (0, 0) fixed point corresponds to the trivial p(r) = 0 = m(r). 


4 . 2 . Power series in terms of the radius 
4.2.1. Linearization about the fixed points 

We suspect that, as in the Newtonian case, the pressure and the mass will have a 
series expansion in the form of an iterated power law. In order to determine the 
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exponents of the power, let’s analyze each fixed point in more detail. Linearizing 
about the (a, a) fixed point we obtain the following system 


where 


Gjlinear 


dC 

d t 


= M + Bx, 


dx 

dt 


= C- x> 


* -2/3 2 , p _ ~2 (l + 5/3 2 ) 

1 + /? 2 ’ (1 + / 3 2 ) 2 


(138) 


(139) 


and G and \ are small perturbations. Therefore, close to the stagnation point, 
solutions are of the form 



A B 
1 -1 


{t - U) 


(m\ 


(140) 


The exponents of the power law are given by the eigenvalues of this matrix, which 
are 


6± 


-1 - 3/3 2 ± v/-7 - 42 /3 2 + (3 A 

2(1 + / 3 2 ) 


(141) 


and the approximate linearized solutions approach the critical point exponentially 
in t — _as exp(i5t) = r s . These are the exact same power exponents found by Harrison 
et. atn Note that if we take the Newtonian limit c —> oo, i.e., a —> 1 and f3 —> 0, 
we recover the eigenvalues for the Newtonian case. 

Let us now analyze the (0, 0) fixed point. Making the substitution x = X an< l 
G = G, where x and G are 1 ° be viewed as small perturbations, we obtain the 
following system 


jGjlinear 




dx 

d t 


(142) 


= G~x- 


Close to the stagnation point, solutions to this system have the form 


(Git) 

\x(t) 


= exp 


0 2 
1 -1 


(t ~ U ) 


(G(U)\ 

\x(ti)J' 


(143) 


This is the second unstable stagnation point, which was mentioned in the previous 
section. This fixed point is a repeller with eigenvalues ( 71 , 72 ) = (2, —1). Note that 
these are the same eigenvalues found for the fixed point (0, 0) in the Newtonian 
analysis. We have now determined the exponents of the critical power-law behaviour 
for both fixed points. 
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4.2.2. Power series about (0,0) in terms of the radius 

Based on Theorem |2] and our general argument for the existence of Frobenius-like 
series let us now assume an expansion for x and C of the form 


OO OO OO OO 

x(r) = E E b mn ?' m7l+ ” 72 ; <(r) = E E a ^n r m7l +" 72 . (144) 

m—0 n —0 m =0 n —0 

We now need to calculate the actual coefficients b n and a n of the series. In order 
to find them, we need to first decouple the system Si. Solving for £ in the d\/dr 
equation in Si and inserting this into the df/dr equation we obtain 



r-i-'X ~ 2y 

dr 


d _ 0 (l + P 2 ) 

dr X 2 1 — 4/3 2 y 


X + r ~r~X 
d r 


X (1 +p 2 ) +/3 2 r^:X 


.(145) 


It is now convenient to use the dimensionless variable z = ry/^Ppo/c 2 . Rearranging 
the above, we obtain 


dz d z 


X 


„ d d d 1 . 

(146) 

2 (1 + 6/3 2 + (3 4 ) + (1 + 5/3 2 + 2/3 4 ) zj^ - 2p 2 z^-z-^- 

+20 1 (3 + /3 2 ) (^x 

Note here that x(°°) = a an d that if we take /3 —> 0 we recover the Newtonian 
result x(oo) = 1. We now have a second order decoupled differential equation for x- 
Before solving for the coefficients of the series, we remind ourselves that for this case 
the 1/r behaviour is unphysical. Hence, let us suppress all occurrences of 72 = —1. 


and restrict our analysis to 71 = 2. Instead of equation (144) we now have 


x(z) = E bn z2n 

n=1 
00 

c( z ) = E® n z2n 

n =1 



(147) 

(148) 


Plugging this ansatz into our decoupled differential equation, we obtain the follow¬ 
ing recursion relation for the coefficients b n 

u _ (l + 6 (3 2 + /3 4 ) + (l + 5 /3 2 + 2/3 4 ) n - d/3 2 n 2 ^ u u 

bn — 7T7o —T\ / . ViOn—i 

(n — 1) (2n + 1) ^ 


4 13 2 (3 + /3 2 ) 
{n — 1) (2 n + 1) 


n— 1 

E i{n — i)bi bn-i- 
z— 1 


(149) 


With this recursion relation, we have an expression for xi z ) equation (147) in the 
form of a power series close to the fixed point ( 0 , 0 ). 
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A similar power series will exist for £. Using the dm/dr equation from Si, we 
observe that (as in the Newtonian case) 

a n = (2n + 1) b n . (150) 

In order to completely specify the power series, we need to have boundary con¬ 
ditions. The boundary conditions to begin the recursion relation are b o = 0 and 
b, = 1/6, which are the same conditions used previously in the Newtonian analysis 
of the system. Choosing b\ = 1/6 is equivalent to the specific choice of z made 
above. Hence, inserting these conditions into the a n equation we obtain cio = 0 
and a\ = 1/2. We have now a complete expression for x and C close to the (0,0) 
stagnation point given by equation ( |147|) and ( |l48| ). 

The radius of convergence R((3) of the \ series will be the same as that for 
the f series, and is now a function of the parameter (3. We performed numerical 
calculations and discovered that there is a finite radius of convergence for almost 
all values of (3. In order to better understand the convergence of the series, let us 
define two regions: Region I, ranging from 0 < (3 < 0.0722 and Region II, defined 
from 0.0722 < f3 < 1. 

The radius of convergence is well-behaved as a function of (3 in region I, allowing 
us to recover the Newtonian upper bound at (3 = 0 (R( 0) ~ 3.2.) Our results 
also showed that the radius of convergence strictly increases monotonically until it 
reaches (3 ~ 0.0722 at which i?(0.0722) « 3.6736. 

In region II, the global behaviour of (3 is given by a power law (approximately 
R{/3) oc /3 _1 ), decreasing from R ss 3.6736 at /3 ~ 0.0722 to R ~ 1.1218 at (3 = 1. 
The local behaviour as a function of (3 in region II, however, is more complex, since 
there are values of [3 for which the radius of convergence drops to zero. The initial 
value at which our series diverges is (3q = 0.0724, and after that it will diverge at 
(3 = (3 o + to, where n is an integer (n = 0,1, 2, 3,...) and r ss 0.00026. It seems 
possible to approximate the behaviour of the radius of convergence in region II with 
a damped oscillating function, e.g. R(f3) oc cos(27 t(3/t + 4>a)/(3- 


4.2.3. Power series about (a, a) in terms of the radius 

In order to determine the coefficients b mn and a mn for the power series about the 
fixed point (a, a) we need to make a shift to \ — X ~ a i n our decoupled differential 
equation (|145| ). By doing so, equation (145) can be simplified and put in operator 
form to obtain 


cl d 


(1 - 2a [1 + 5 f3 2 + 2/3 4 ]) z- - (l - 4 a/3 2 ) z-z- - 2 


dz d z 


X 


d d 


2 (1 + 6f3 2 + /3 4 ) + (1 + 5 f3 2 + 2/3 4 ) 2 — - 2(3 2 z—z— 


dz dz 


X 


+2/3 2 (3 + /3 2 ) (4x) 2 - 


(151) 
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Again, we have a second order decoupled differential equation for Note that for 
our series expansion we have the expression given in equation (144), where we must 
take Ai ,2 —> + 2 - Plugging this ansatz into our decoupled differential equation, we 
obtain a recursion relation for the coefficients b mn 


2 (1 + 6 ■0 2 + p 2 ) + (1 + 5/3 2 + 2/3 4 ) A - 2+ A 2 ^ ^ ^ 

“ ~ (1 - 4/3 2 a) A 2 - (1 - 2a [1 + 5/3 2 + 2/3 4 ]) A + 2 EE b H b rn-i-,n-j 
v v L 0 2=0 

2+ (3 + /3 2 ) 

_ (1 - 4/3 2 a) A 2 - (1 - 2a [1 + 5/3 2 + 2/? 4 ]) A + 2 

m n 

*EE ij (to - j) ( n-i)b ij bm—i',n—j 5 (152) 

j =0 i =0 


where A = mdi + 7182 - With this recursion relation, we again have a complete 
expression for \ iu the form of a power series close to the fixed point (a, a). ( will 
also have a power series expansion, with coefficients given (as in the Newtonian 
case) by 


ttmn — (A ~\~ 1) b mn . 


(153) 


Once more, we need boundary conditions to begin the recursion relation. We know 
analytically that 600 = 0 and a 00 = 0. If we knew &01 = ^*0 [° r equivalently 
aoi = al 0 = (1 + c+froi] then we could recursively calculate the entire series. As in 
the Newtonian case, there is no analytic theory for these coefficients. 

Though we do not have an analytic way of specifying aoi or 601 > we do have 
the powerful result that C and \ are [close to the fixed point (a, a)] given by power 
series of the form 


and 


OO OO 

x (r) = a + E E 

m=0n—0 


/ 47rp 0 r 

l c 2 


m 8 \-\-n 82 


OO OO 

CM =q + EE amn 

m=0 n=0 


( , \ mS 1 +nS 2 

47rp 0 r j 


(154) 


(155) 


In terms of the dimensional physical variables p(r ) and m(r) this can explicitly 
be put in Frobenius-like form 


and 


i(r) = 2c 2 r 


p(r) = 


OOOO / , 


771=0 71=0 


771 <5l +71(52 


47T7-2 


OO OO 

a+EE ^77171 

771=0 71 = 0 


/ 47rpo^ 


771 ^ 1 + 71(52 


(156) 


(157) 
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4.3. Power series in terms of the pressure 

Though we will not go into any of the details, it is clear from the existence of the 
power series p(r), that the series can be “reverted” to provide a power series for 
r(p). The key points are that p(r) is monotonic, so the inverse r{p) exists. Then 
from the power series above we deduce the existence of a reverted Frobenius-like 
power series 



(158) 


From this we deduce that for a cutoff equation of state (linear in the core, constant 
density in the envelope) even in the relativistic case the core mass and core radius 
will exhibit damped oscillations as a function of central pressure po- 



Similarly the total mass and total radius will exhibit damped oscillations. 

It should also be clear that the details of the equation of state at low density 
do not matter. As long as the equation of state at high densities is linear, then the 
analysis above implies that as the central pressure goes to infinity the core radius 
and core mass will undergo damped oscillations as a function of central pressure. 
Once the core undergoes these oscillations, the envelope serves only to communicate 
this internal behaviour out to the stellar surface. Thus the total mass and total 
radius will similarly undergo damped oscillations. These damped oscillations will 
involve an infinite power series in terms of p j^, which will contain the first-order 
term discussed in Harrison et ala , plus an infinite collection of higher-order terms. 

5. Conclusions 

We have presented a unified formalism for calculating the solution to scale invariant 
theories in the form of power laws, and deviations from pure power-law behaviour. 
We have proved that some of the solutions to a scale-invariant differential system, 
whose associated autonomous equation possesses fixed points, will exhibit pure 
power-law behaviour. We have also proved that after linearizing about these fixed 
points, the solution to the scale-invariant linearized system can also be written in 
the form of an independent collection of power laws. Generically, these collections 
of power-law solutions can be combined to develop a Frobenius-like series in a 
neighborhood surrounding the fixed point. 

We applied these ideas to several problems taken from astrophysics, such as 
static isothermal stars both in the Newtonian and TOV frameworks. We showed, 
both numerically and analytically, that in all the static cases the solution to the 
scale-invariant differential system associated with each example exhibits power-law 
behaviour as the star approaches the point of collapse. We also showed analytically 
that the damped oscillations in the pressure, mass and radius of a star, as the 
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central density goes to infinity, can be represented as a Frobenius-like power series, 
which is easily obtained by applying our formalism. 

More generally, it would be advantageous to consider truly dynamical systems 
involving PDEs in both space and time. Work along these lines is progressing. 
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Appendix A. Autonomous equation without fixed points 

Caution must be exercised since not all scale invariant equations possess fixed 
points. A scale invariant differential equation, whose associated autonomous equa¬ 
tion does not have any fixed points, will not possess power-law solutions. As an 
example, let’s analyze what happens to the following relatively simple differential 
equation: 

r-i/ / ^ _ d y py{x) c „p-i n nil 

F(x, y(x)) = - f x p =0, (A.l) 

da: x 

where / is a non-zero constant. This is easily verified to be a scale-invariant equation 
and the associated equidimensional-in-a; equation is 

F(x,w{x)) = x p =0. (A.2) 

\ ax x J 

This factorizes to give 

m nm = - /) = o. (a.s) 

This last equation is autonomous and it does not possess a fixed point since we are 
assuming that / ^ 0. The general solution for this differential equation is easily 
verified to be 

y(x) = x p (/ In x + A). (A.4) 

Note that the solution to our example will be a power law if and only if / = 0, 
regardless of the fact that our example is explicitly scale invariant for all values of 
/. In conclusion, scale invariance by itself does not guarantee that the solution to a 
scale invariant differential equation will exhibit power law behaviour. More compli¬ 
cated examples of this behaviour can easily be constructed, but this simple equation 
is already enough to make the point that scale invariance does not necessarily lead 
to power-law solutions. 

However, if we impose both scale invariance and the constraint that the associ¬ 
ated autonomous equation has fixed points, then the solution will certainly exhibit 
power law behaviour. 
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Appendix B. Limit cycles and discrete self-similarity 

In the theory of autonomous differential equations the existence of a fixed point is 
one of the simplest things one could wish for; the next most complicated structure 
one might encounter is a limit cycle. 

Suppose we start with a scale invariant differential equation and that the asso¬ 
ciated autonomous equation has a limit cycle rather than a fixed point. This means 
there is a special solution z*(t) with a period T such that 

z*(t + nT) = z*(t) VneZ. (B.l) 

In terms of the equidimensional variables this implies 

w*(exp(nT)x) = w*(x). (B.2) 

That is, a limit cycle in the variable t becomes a discrete self-similarity in the 
variable x. Similarly in terms of the original variable y(x) we obtain discrete power- 


law behaviour 

y(x) = x p w*(x); Vn € Z : y(exp(nT)x) = exp (pnT) y(x). (B.3) 

This is not the classic power-law behaviour discussed previously, but is, in a sense, 
the next least complicated thing. As an example, let’s consider the second-order 
differential equation 

F(x, y{x)) = x 2 y"{x) — (2 p — 1 )xy'(x) + p 2 y + O 2 y(x) = 0 (B.4) 

This is easily verified to be scale invariant with index p. The associated equidimensional- 
in-a; equation is 

F(x, w(x)) = x 2 w"{x) + xw'(x) + il 2 w(x) = 0. (B.5) 

The associated autonomous equation is 

F(z(t)) = z{t) + fl 2 z(t) = 0. (B.6) 

Now there is a trivial fixed point at z = 0 but the general solution to the autonomous 
equation is 

z(t) = Acos(Clt) + Bsin(fit). (B.7) 

In terms of the equidimensional variables 

w{x) = 9? [{A + iB)x in ] . (B.8) 

Finally, in terms of y{x) 

y(x) = x p $ [(A + iB)x m ] . (B.9) 


This now is an explicit example of a discrete power law. The solution is not scale 
invariant under arbitrary rescalings x —> a x, but it is invariant under the specific 
rescaling x —> x exp(27r/fl). 

Note the implication: Discrete self-similarity and discrete power laws are not 
as peculiar as one might at first imagine; instead they are simply the next most 
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complicated thing (after simple power-law behaviour) that can happen in scale 
invariant systems. (Viewed i_n_this light the occurrence of discrete self-similarity in 
Choptuik’s critical solution Lr0 should not at all be considered surprising.) 


Appendix C. Isothermal Lane Emden equation. 

The traditional [isothermal] Lane-Emden equation is 




> A(p(z) 

dz 


= exp {-(j){z)). 


m = 0 ; g 


= 0 . 


(C.l) 


(C.2) 


This second-oder DE is completely equivalent to the first-order system we derived 
for a Newtonian relativistic star. Indeed 


</>(z) = - In [ 2 C(^)/- 2 ] , 


(C.3) 


and the series we previously determined for £(z) in the Newtonian case then carry 
over into series for the Lane-Emden function (j>(z). It is then straightforward to 
compute 


<j)(z) = -z 2 -—z 4 + — 
^ ’ 6 120 1890 


z 6 - 


61 


629 


1632960 


224532000 


z w + 0(z 12 ). (C.4) 


You can also derive this in a more direct way by putting a trial power series directly 
into the Lane-Emden DE and equating coefficients. 

Perhaps more surprising is the implication that at large z the Lane-Emden 
function should have an expansion of the form 


<Kz) = ln[z 2 /2] z^ mSM \ (C.5) 

m—0 n=0 


[Recall Si t 2 = (—1 ± iy/7)/2.] The coefficients q mn are in turn determined by the 
recursion relations for a mn and b rnn previously discussed. An alternative more direct 
route is to write 


<t>( x ) = ln[z 2 /2] + q(z), 

and substitute into the Lane-Emden equation yielding 

{ Z ~^z + q ^ = 2 { exp [ _9 ( z )] “ ' ( C - 6 ) 

Linearizing around the obvious solution q(z) = 0 recovers the critical indices = 
(—1 ± i\] 7)/2, while inserting a trial series consisting of integer powers of these 
critical exponents will determine all the higher coefficients q mn in terms of <701 = 9io ■ 
Of course, as is usual, there is no analytic theory for this first coefficient since it 
corresponds to finding the unique boundary condition at infinity that leads to a 
regular star at the center — goi = 9io has to be determined numerically from an 
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outward integration of the Lane-Emden equation starting from the origin where 
we do know the physical boundary conditions. 

When it comes to writing down explicit recursion relations for the coefficients, 
it is more convenient to work with the compactness x( z ) i n terms of which the 
[isothermal] Lane-Emden equation is equivalent to equation ( f78| ) and explicit re¬ 
cursion relations are given in equations (|8(]) and (jS9|). 
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